SoSe2021

Folienübersicht

(Nicht-parametrische) Mehrstichprobentests

Mittelwerte von ≥ 3 Stichproben vergleichen

Das Problem der multiplen Vergleiche

Beispiel

  • Auswirkung der Meerestiefe auf das bakterielle Wachstum
  • Probenahme in 5m (Kontrolle), 10m, 20m, 40m, 80m, 160m

→ Warum nicht jede Tiefe mit der Kontrolle mittels t-Test vergleichen?

‘Familywise Error Rate’ (FWER)

  • Sehr ineffizient!!! 5 Tests sind notwendig
  • Je öfter ein Test mehrfach auf eine Hypothese angewendet wird, desto höher ist das Risiko eines stat. Fehlers 1. Art (\(\alpha\); Ablehnen einer korrekten H0).
  • Die FWER ist dabei definiert durch die Wahrscheinlichkeit, mindestens eine der Testhypothesen fälschlicherweise abzulehnen:
    • \(\alpha_{FWER} = 1-(1-\alpha)^c \Rightarrow 1-(1-0.05)^5 = 0.23\) (c = Anzahl der Vergleiche)
    • => Der gemeinsame Fehler ist fast fünfmal höher als der nominell gewählte. Die Wahrscheinlichkeit, dass mindestens ein Vergleich zu einem stat. Fehlers 1. Art führt, beträgt 23%.

Lösungen zur Kontrolle der FWER

Paarweise Tests

  • Bonferroni-Korrektur
    • Kann immer angewendet werden, jedoch recht konservativ → Verlust an Teststärke
    • Leicht zu berechnen: \(\alpha_{adj} = \frac{\text{Signifikanzniveau}~\alpha}{\text{Anzahl der durchgeführten Tests}} \Rightarrow \frac{0.05}{5} = 0.01\)
    • Alternative Variante: Bonferroni-Holm-Verfahren (Sequentielle Bonferroni-Methode mit höherer Teststärke)
    • Funktion p.adjust() wandelt unkorrigierte p-Werte in korrigierte mittels Bonferroni-, Bonferroni-Holm- und weiteren Korrekturen um.
  • Spezifische Tests für multiple Mittelwertvergleiche (sog. Post-hoc-Tests):
    • Tukey’s, SNK, Ryan’s etc.
    • passen Signifikanzniveaus auf unterschiedliche Weise an

Gruppenweise: Mehrstichprobentest

  • Parametrisch:
    • Varianzanalyse (ANOVA)
  • Nicht-parametrisch:
    • Kruskal-Wallis-Test (H-Test)
    • Friedman-Test (Friedman-Rangsummentest)

Verteilungsfreie Tests

Kruskual Wallis Test (H Test)

  • Logische Erweiterung des Mann Whitney U-Tests, Messung kann also auch auf einer Ordinalskala erfolgen
  • Rang-basiertes Gegenstück zur 1-faktoriellen ANOVA → allerdings nur 95% der Teststärke
  • Stichprobengröße der Gruppen kann unterschiedlich sein
  • Annahme: Stichproben sollten annähernd die gleiche Verteilung haben

Friedman Test

  • Rang-basiertes Gegenstück zur 1-faktoriellen ANOVA mit wiederholten Messungen
  • Messung kann auch auf einer Ordinalskala erfolgen
  • Anwendbar auf vollständige Blockdesigns
  • Teststärke oder Effizienz: abhängig vom k/n-Verhältnis

Teststatistik H

\[H=\frac{12}{N*(N+1)}*\left( \sum\frac{R_i^2}{n_i} \right)-3*(N+1)\]

N = Gesamter Stichprobenumfang, \(n_i\) = Anzahl Messung pro Gruppe \(i\), \(R_i\) = Summe der Ränge pro Gruppe \(i\)

Teststatistik \(F_r\)

\[F_r=\left( \frac{12}{n*k*(k+1)}*\left( \sum R_i^2\right)\right)-3*n*(k+1)\]

n = Anzahl der Blöcke, k = Anzahl der Gruppen, \(R_i\) = Summe der Ränge pro Gruppe \(i\)

Kruskual Wallis Test

Einzelschritte

  1. Die Rohdaten aus den k Stichproben/Gruppen werden zu N Elementen zusammengefasst und dann vom niedrigsten zum höchsten Rang geordnet.
  2. Gleiche Werte werden mit dem Mittelwert der jeweiligen Ränge versehen.
  3. Diese Ränge werden dann wieder in die k Gruppen einsortiert.
  4. Bildung der Rangsummen pro Stichprobe i=1 bis k
  5. Berechnung der Teststatistik H
    • Wenn mehr als 25 % den gleichen Rang haben, muss die H-Statistik korrigiert werden
  6. Prüfverteilung
    • wenn \(n_i\) < 6 und k = 3 wird die H-Test-Tabelle verwendet
    • wenn \(n_i\) ≥ 6 und k > 3 ist H näherungsweise \(\chi^2\)-verteilt (df = k-1)

Kruskal-Wallis-Test in R: manuell

Beispiel iris

Rangsummenbildung

df <- iris %>%
  mutate(ranks_sl = rank(Sepal.Length, ties.method = "average")) %>%
  group_by(Species) %>%
  summarise(sum_ranks = sum(ranks_sl))
df
# A tibble: 3 x 2
  Species    sum_ranks
  <fct>          <dbl>
1 setosa         1482 
2 versicolor     4132.
3 virginica      5710.

Teststatistik H und p-Wert

(h_sl <- 12/(150*(150+1))*sum(df$sum_ranks^2/50) - 3*(150+1))
[1] 96.76096
pchisq(h_sl, df = 3-1, lower.tail = FALSE)
[1] 9.741475e-22

Kruskal-Wallis-Test in R: kruskal.test()

Beispiel iris

Vektorsyntax

# x = abhängige Variable, g = Gruppierungsvariable/Faktor
kruskal.test(x = iris$Sepal.Length, g = iris$Species)

Formelsyntax

# kruskal.test(x ~ g, data)
kruskal.test(Sepal.Length ~ Species, data = iris)
    Kruskal-Wallis rank sum test

data:  Sepal.Length by Species
Kruskal-Wallis chi-squared = 96.937, df = 2, p-value < 2.2e-16

Friedman Test in R

friedman.test()

Syntax

friedman.test(y, groups, blocks)
# y = abhängige Variable
# groups = Messwiederholung -> eigentlich zu untersuchender Faktor
# blocks = zu untersuchende Objekte 

# oder Formelsyntax
friedman.test(y ~ groups|blocks, data)

Your turn …

Quiz 1 zum Nachmachen

Kronblattweite in iris

Es soll die Hypothese getestet werden, dass sich die durchschnittliche Kronblattweite der drei Iris Arten unterscheidet.

kruskal.test(Petal.Width ~ Species, data = iris)
    Kruskal-Wallis rank sum test

data:  Petal.Width by Species
Kruskal-Wallis chi-squared = 131.19, df = 2, p-value < 2.2e-16

Varianzanalyse (ANalysis Of VAriance)

Was ist die Varianzanalyse oder ANOVA?

  • In seiner einfachsten Form die Erweiterung der t-Statistik für mehrere Stichproben.
  • Eine weit verbreitete statistische Methode.
  • Testet auf Mittelwertunterschiede, indem die geschätzte Gesamtvariabilität in den Daten in 2 Teile zerlegt wird:
    • Variabilität innerhalb (‘within’)der verschiedenen Gruppen
    • Variabilität zwischen (‘between’) den verschiedenen Gruppen
  • Wenn sich die Varianzen zwischen den Stichproben unterscheiden, müssen sich auch die Mittelwerte zwischen den Stichproben unterscheiden!
  • Mit anderen Worten: je größer der Unterschied in der Varianz zwischen den Stichproben im Vergleich zum Unterschied innerhalb der Stichproben ist, desto wahrscheinlicher ist ein Unterschied im Mittelwert zwischen den Stichproben!
  • Wichtige Annahme: Homogenität der Varianzen innerhalb jeder Gruppe (wichtiger als Normalverteilung!!!)

Variablen - und ANOVA-Typen

  • Variablentypen:
    • Unabhängige Variablen (X) sind meist kategorial und werden “Faktoren” genannt, die verschiedenen Kategorien werden “Faktorstufen” oder “Gruppen” genannt → sollten auch in R explizit als factor definiert werden!
    • Abhängige Variablen (Y) sind diskret oder kontinuierlich
  • ANOVA-Typen (abhängig vom experimentellen Design):
    • 1 Faktor = einfaktorielle Varianzanalyse (1-way ANOVA)
    • 2 Faktoren = zweifaktorielle Varianzanalyse (2-way ANOVA)
    • ≥3 Faktoren = mehrfaktorielle Varianzanalyse (multi-way ANOVA)
    • Wenn >1 Faktor
      • geschachteltes (nested) vs. gekreuztes (crossed) Design
      • fixed vs. random vs. mixed effect model
      • mit/ohne Messwiederholung, Blockdesign

Das Regressionsmodell

Das ANOVA Modell

Das Prinzip der ANOVA

1

Das Prinzip der ANOVA

1

Varianzzerlegung

Ein Anwendungsbeispiel: Lachszucht

Gewichtsunterschiede von Atlantischem Lachs, der in 4 verschiedenen Typen von Netzkäfigen gezüchtet wurde.

  • Einfluss der Käfigtypen auf das mittlere Gewicht des Atlantischen Lachses (Salmo salar)
  • Antwortvariable Y: Körpergewicht [in kg]
  • Prädiktorvariable X: Käfigtyp
    • kategorial mit 4 Gruppen: A,B,C,D
  • Replikate sind Individuen pro Käfig i (\(n_i\) = 24)

Erkärte vs. nicht erklärte Varianz

Replikat Typ A Typ B Typ C Typ D
1 2.0 2.7 2.9 2.2
2 2.7 2.0 1.9 2.0
3 2.2 3.8 2.7 3.1
… … … … …
24 3.1 1.5 2.5 2.7
Mittelwert 2.5 2.5 2.5 2.5

Gesamtmittelwert: 2.5

  • Behandlungen (hier = Käfigtypen) erklären nichts
  • SSGruppen = 0

Erkärte vs. nicht erklärte Varianz

Replikat Typ A Typ B Typ C Typ D
1 5.0 4.1 2.9 0.8
2 5.0 4.1 2.9 0.8
3 5.0 4.1 2.9 0.8
… … … … …
24 5.0 4.1 2.9 0.8
Mittelwert 5.0 4.1 2.9 0.8

Gesamtmittelwert: 2.5

  • Behandlungen (hier = Käfigtypen) erklären alles
  • SSResidual = 0

Typische ANOVA Tabelle

Sources of variation SS df MS F p
Groups \(n \sum(\bar{y_i}-\bar{y})^2\) \(p-1\) \(\frac{SS_{Groups}}{(p-1)}\) \(\frac{MS_{Groups}}{MS_{Residuals}}\)
Residuals (or within) \(\sum\sum(y_{ij}-\bar{y_i})^2\) \(p(n-1)\) \(\frac{SS_{Groups}}{(p-1)}\)
Groups \(n \sum(y_{ij}-\bar{y})^2\) \(pn-1\)

H0:

  • Kein Unterschied zwischen den Gruppenmitteln \(\mu_1 = \mu_2 = ... = \mu_i = \mu\)
  • Denn Gruppeneffekte sind definiert als \(\alpha_i = \mu_i - \mu\),
    • Kann auch ausgedrückt werden als: “Es gibt keinen Effekt bestimmter Behandlungsstufen (treatment effects), die Unterschiede sind nur zufällig”: \(\alpha_1 = \alpha_2 = \alpha_3 = .. = 0\)

Sowohl \(MS_{Groups}\) als auch \(MS_{Residuals}\) schätzen \(\sigma^2\) → Quotient = 1

  • F-Ratio folgt F-Verteilung mit \(df_{Groups}\) und \(df_{Residuals}\)

Beispiel: Gewichtsunterschiede bei Atlantischem Lachs

Sources of variation SS df MS F p
Groups 52.10 4-1 171724 17.367 <2e-16
Residuals (or within) 26.02 4*(24-1) 0.283
Groups 78.12 4*24-1
  • Es gibt einen signifikanten Effekt von mindestens einem Käfigtyp auf das mittlere Gewicht des Atlantischen Lachses. Mit anderen Worten: Die mittleren Gewichte von Atlantischem Lachs, der in verschiedenen Käfigen gezüchtet wurde, unterscheiden sich mindestens bei einem Käfigtyp.
  • ABER: wir wissen noch nicht, welche Gruppe sich unterscheidet!

ANOVA in R

aov()

aov(formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, ...)

Die volle ANOVA Tabelle erhältst Du aber nur mit der summary() Funktion, die auf das aov Objekt angewendet wird:

model <- aov(formula, data)
summary(model)

ANOVA in R

Lachsbeispiel

mod <- aov(formula = weight ~ cages, data = salmon)
mod
Call:
   aov(formula = weight ~ cages, data = salmon)

Terms:
                   cages Residuals
Sum of Squares  53.91378  18.51342
Deg. of Freedom        3        92

Residual standard error: 0.4485898
Estimated effects may be unbalanced
summary(mod)
            Df Sum Sq Mean Sq F value Pr(>F)    
cages        3  53.91  17.971   89.31 <2e-16 ***
 [ erreichte getOption("max.print") --  Zeile 1 ausgelassen ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Überprüfung der Modellannahme posteriori

op <- par(mfrow = c(2,2))
plot(mod)

par(op)

Effektgröße \(\eta^2\) bei der 1-faktoriellen ANOVA

plot.design()

Die Berechnung der Effektgröße für Designs zwischen Gruppen ist viel einfacher als für Designs innerhalb von Gruppen:

\[\eta^2 = \frac{SS_{Groups}}{SS_{Total}}\]

mod_sum <- summary(mod)
# str(mod_sum)
eta_sq <- mod_sum[[1]]$`Sum Sq`[1] / (mod_sum[[1]]$`Sum Sq`[1] + mod_sum[[1]]$`Sum Sq`[2] )
eta_sq
[1] 0.7443858

Effektgröße grafisch darstellen

plot.design()

plot.design(weight ~ cages, data = salmon)

Fragen?